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At low temperatures and voltages tunneling transport through an interface between a supercon- 
ductor and hopping insulator is dominated by coherent two-electron tunneling between the Cooper- 
pair condensate and pairs of localized states, see Kozub et al, Phys. Rev. Lett. 96, 107004 (2006). 
By detailed analysis of such transport we show that the interface resistance is extremely sensitive to 
the properties of the tunneling barriers, as well as to asymptotic behavior of the localized states. In 
particular, dramatic cancellation takes place for hydrogen-like impurities and ideal barrier. However, 
some disorder can lift the cancellations restoring the interface transport. We also study non-Ohmic 
behavior of the interface resistor and show that it is sensitive to the Coulomb correlation of the 
occupation probabilities of the involved localized states. It is expected that non-Ohmic contribution 
to I — V-curve will experience pronounced mesoscopic (fingerprint) fluctuations. 

PACS numbers: 72.20.Ee, 74.45. +c, 74.45. +r 



I. INTRODUCTION 

In this paper, we address the charge transfer through 
the interface between a superconductor (SC) and a hop- 
ping insulator (HI), i.e., a system where transport oc- 
curs via hops between localized (non-propagating) single 
particle states. There are many experimental situations 
in which the hopping insulator is coupled to a measur- 
ing circuit via superconducting leads, see, e.g., Ref. Q. 
However, it has long been known that transport of sin- 
gle electrons into or out of a superconductor is expo- 
nentially suppressed at low temperatures as e~ A / T , due 
to the energy gap, A, of the superconductor^ Conse- 
quently, single-electron tunneling cannot be responsible 
for the charge transfer between a superconductor and any 
normal conductor. 

The problem of charge transfer between a SC and a HI 
was first addressed in Ref. 0, where it was shown that at 
low temperatures the transport is governed by the time- 
reversal reflection, where electrons tunnel into supercon- 
ductor from the localized states in the hopping insulator 
located near the interface. This process is similar to so- 
called crossed Andreev charge transfer discussed previ- 
ously in connection with various mesoscopic systems^^ 
Electrons hopping from the superconductor to impuri- 
ties near the surface of the insulator and back again have 
been proposed as a source of qubit decoherence for some 
systems £ 

In this paper we extend the analysis of Ref. in two 
directions. Firstly, we consider the influence of the prop- 
erties of the tunneling barrier on the charge transfer, 
which turns out to be surprisingly sensitive to the bar- 
rier roughness. Namely, we found that interference ef- 
fects in tunneling can lead to a significant increase in 
the interface resistance due to fine cancellations of the 
contributions to the two-particle tunneling probability. 



The roughness of the barrier suppresses these effects and 
in this way influences the interface resistance. Secondly, 
we consider non-Ohmic transport through the interface. 
We will show that the interface contribution to the resis- 
tance can be strongly nonlinear, and that the nonlinear 
behavior is essentially related to the Coulomb correlation 
of the occupation numbers of the localized states adja- 
cent to the interface. One can expect a rich pattern of 
reproducible (fingerprint) fluctuations in the I — V-cuive 
due to pronounced non-Ohmic contributions of individ- 
ual pairs. Thus, combined studies of the linear interface 
resistance and average nonlinear I —V curve and its fluc- 
tuations may tell a lot about barrier details and about 
the formation of the depletion zone near the barrier. 

The paper is organized as follows. In sec. [H] we set 
the stage presenting the model of Ref. for the coherent 
charge transport. Detailed calculations for the Ohmic 
case are presented in Sec. Ill where we show how the 
interference-based cancellations occur. In Sec. llVl we dis- 
cuss several generalizations of the model, including the 
effects of modification of the tails of localized wave func- 
tions and barrier details. In Sec.0 starting from calcu- 
lation of the non-Ohmic contributions of the individual 
pairs, we demonstrate how account of the Coulomb cor- 
relations leads to non-Ohmic behavior of the interface 
conductance. 



II. MODEL 

As shown in Ref. 0, the contact resistance can be gov- 
erned either by the interface tunneling barrier, or by de- 
formation of the hopping cluster in the HI in the vicinity 
of the interface. Here, for brevity, we will focus on the 
situation where the interface resistance is dominated by 
the barrier. 
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We start with the case of linear conductance where 
it is natural to use the Kubo linear response theoryi 
According to this theory, the conductance, G, is expressed 
through the susceptibility, 



jt(i),/(0) 



dt 



(1) 



as Q = lim^^o LU ~ 1 Im^(w). Here I(t) is the current 
operator and we will use units where h = 1. 

Let a superconductor and a hopping insulator occupy 
adjacent 3D semi-spaces separated by a tunneling bar- 
rier (TB). The presence of the barrier simplifies calcula- 
tions which will be made in the lowest non- vanishing ap- 
proximation in the tunneling amplitude To. This models 
the Schottky barrier usually present at a semiconductor- 
metal interface. Then the current operator is defined as£ 

where r is the coordinate on the superconductor side of 
the interface plane, r' is the coordinate on the semicon- 
ductor side, a + (r, t) and &(r, t) are creation and annihila- 
tion operators in the semiconductor and superconductor, 
respectively, d is the electron localization length under 
the barrier. 

The Matsubara formalism, see e.g. Ref. 0, allows cal- 
culation of the susceptibility by analytical continuation 
of the so-called Matsubara susceptibility defined as 



Xm(O) 



T t /(t)/(0) 



3 lUT dT . 



Here T T means ordering in "imaginary time" , r, (3 = f /T, 
temperature T is measured in units of energy. The inte- 
gration over the imaginary time actually means the av- 
erage over a grand canonical ensemble with temperature 
T and chemical potential fj,. 

The operators P and I are time-dependent interaction 
picture operators. Changing to Schrodinger type opera- 
tors we write 

XM {Q) = J (T T I(r)I(0)eJ HTdT ) e mr dr 

where Ht is the tunneling Hamiltonian given by the ex- 
pression 



H t (t) = J d 2 rdVT(r,r')[a t (r ! r)6(r',r)+h. 1 



(2) 



where the integration is performed along the interface. 
Here T(r, r') is the tunneling amplitude which in general 
is dependent on the coordinates both for entry to and 
exit from the barrier. 
Let us first assume that 



as it was done in Ref. y. Then 

H t (t) =T J d 2 r[aHr,T)b(r,T)+h.c.} . (4) 

Because single electron transitions are forbidden by the 
energy gap, we have to expand the expression for x to the 
second order in Ht- Of the many possible contractions, 
we are only interested in two-electron transitions in both 
directions. We get a total of 12 different contractions. 
Of these, half will be only the hermitian conjugate of 
the other half, and a symmetry consideration reduces the 
number of relevant contractions to 3. They are 

[a\.b T - a T bl][alb - a bl][a\b 1 +aib\][alb 2 +a2bl], (A) 



[ atb T - a r 6t][4bo-a 6j ) ][a t 1 & 1 + ai b\] [a|&2+a 2 j >|], (Bl) 

[ a f T b T - a T 4][a| ) bo-a & ( ! ] ][a|& 1 +a 1 6|][a^ 2 +j2^ ]- ( m ) 

One can show that the first one (A) is small at eV <C A, 
i.e., when the single-electron transport is suppressed, 
while the two others, Bl and B2, give equal contribu- 
tions. We will therefore follow only Bl through the fur- 
ther analysis. To perform calculations we introduce the 
Green's functions in a usual way 1 ' 

(T r 6(r,r)fe(r',r')) = F{x-x'), 
(T T tf(r,T)tf(r',T')) = F\x-x'), 
(T T a(r,T)a^(r',T')) = G(x,x'). 

Here Xi = (rj,r,), F(x — x') is the anomalous Green 
function in the superconductor, and G(x, x') is the Green 
function of the insulator. We get 



T t I{t)I{0)J = e z T^ J d z rd% dxi dx 2 

xF(x - X!)F^(x 2 - x )G(x ,x)G(x 2 ,x 1 ) . (5) 

Here dx = (fir dr. We then take the discrete Fourier 
transforms of the Green's functions to write them as func- 
tions of f2, uji rather than r, Tj . In this case we will have to 
require that fl and u> take discrete values of Q m = 2i:mT 
and ui n — (2n + \)nT. Performing the integrations over 
t, n and t 2 , we get (5-functions for the relations between 
the different discrete frequencies, giving 



x(n m ) = 2Te 2 T*J2 Jll d2 



xF(r - n^njf^ro - r 2 , -w„) 

xG(r ,r, uj n - Q m )G(r 2 , ri, -w n ) . (6) 

The Feynman diagram corresponding to this expression 
is shown in Fig. ^ 

We assume the localized states to have hydrogen-like 
wave functions, centered on impurities in positions r s , 
with energies e s and localization length a, 



T(r,r')=T J(r-r') 



(3) 



tf,(r) = (^a 3 )- 1 / 2 e-l r - r =l/ a 



(7) 
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and the insulator Green's function is 

*:(r)* 8 (r') 



G(r,r>„)=]T 



(8) 



We now define r s and 17 as coordinates of the impuri- 
ties contributing to G(r 2 ,ri) and G(r ,r), respectively. 
Then the Matsubara susceptibility can be expressed as 



For the anomalous Green's function we use 
d 3 p A 



F(R,. n ) = f* 



— ip-R,/fr 



(2nh)3 A 2 + Q 



ng m A sin{Rk F ) R V**+»l 



2^ 



Rki 



(9) 



Here £ p = (p 2 — p F )/2m, £ is the coherence length in a 
superconductor, g m = mp F /ir 2 is the density of states in 
a metal. 



III. CALCULATIONS 

So far we just followed Ref. 0, but we will now demon- 
strate how the oscillations of the anomalous Green's 
functions, oc sinRkp, lead to a significant decrease of 
the result comparing to the simple estimates presented 
there. It turns out that these oscillations give rise to pro- 
nounced cancellations in the susceptibility for the case of 
a hydrogen-like impurity state. 




FIG. 1: Diagram relevant for the Matsubara susceptibility. 
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E 



A 2 + lo 2 (— iu} n - e s )(iw n - iVt m - ej) 



^^ shrfcHr-nl 
k F \r - ri 



|r-r,| (Ir-nl)VA 5 



<A 



(10) 
(11) 



It is safe to assume that |r — ri | and |ro — r%\ are of the 
order the distance p s i between the impurities projected 
onto the interface, but with variations of the order a, 
where kpa 1. The localization length a can be esti- 
mated as the Bohr radius do = AttkH 2 /m*e 2 . Assuming 
m* rs 0.1m e where m e is the mass of a free electron, and 
K ~ 10, we get k F a sn 100. Since the superconductor 
localization length is much greater, ( > a, we can safely 
replace |r — ri| — > p s i in Eq. (|11|) in all places except 
the sinfcp|r — ri|, which oscillates rapidly. Choosing the 
coordinates as shown in Fig. [21 we obtain 



k F p s i 

x J d 2 rd 2 n sinfc F |r-ri|e" (|ri - rs|+|r - ri|)/a . 

As r and are located in the interface plane, it is 
natural to choose the origins for r and to be at the 
projections of r s and 17, respectively, into the interface 




FIG. 2: Coordinates for the spatial integration. 

plane of the. Then d\% = Vi dr,; d9i, 

|r - rj| = yjrl+zf, |n - r s | = ^ r 2 + z 2 s , |r x - r 2 | 
= \J {psi + rsm8 — r\ sin^i) 2 + (r cos9 — r\ cos#i) 2 

where z s and zi are the distances from the impurities to 
the interface. 

The minimal distance between the impurities tak- 
ing part in the coherent transport, p"" n , is limited by 
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Coulomb correlation. If a/ pff 1 « (kpa) -1 one can ne- 
glect the items containing cos 9i in the above expression. 
Then |r — i"i | r* p s i + r sin 9 — n sin 9\ , and the integrals 
over r and ri can be calculated exactly using the formula 



d9 sin(k F r sin# + <p 

-7T 

Ime'* f d9e tkFrsin 



= Im / d9e lkFrsin9+l ' 1 ' 



27rsin(0) Jo(fcr) . 



After integrating over both angles we get 
(2tt) 2 sin kppsi ^ 



hi = 
I r (z) = 



kFPsl 

rdrJo(k F r)e 



3 I r (z s )I r ( Zl ) 



(12) 



The integral I r will be referred to extensively in later sec- 
tions, as most of the modifications we will discuss change 
this integral only, leaving all other calculations unaltered. 
We can then use the following identity^ 



xdxe- pVx " +z \j (cx) 



= p(p 2 + c 2 )- 3 / 2 (l + z^T^) e"*V^ 
[Re(p) > |Im(c)|; Re(z) > 0] . 

Combining all integrals, this gives: 

(27r) 2 a 4 sinfc F p s! T{zi)T{z s ) 



hi = 

Hz) = 



kFPsl 

z 



(l + £; 2 a 2 ) 3 



1 + k 2 F a 2 | 



Using the assumption that kpa. 3> 1 the function T{z) 
simplifies to 



T[z) w (1 + k F z)e 



The essential observation here is that the expression for 
I s i contains a factor e ~^ Zl+Za ' >kF , which again should be 
squared for the final result. This conclusion contradicts 
Ref. 0, where the factor e~ < - z,+Zs ^ a was predicted. The 
strong decay of hi as a function of zi and z s means that 
only pairs very close to the interface can contribute to 
the Cooper pair transfer. This means that the theory as 
presented above and in Ref. proves its own inadequacy, 
since the assumption that the wave functions of the local- 
ized states is bulk hydrogen-like ones requires z\ and z s 
to be at least of the order of a. For closer impurities, the 
wave function is modified by the vicinity of the surface, 
and the result becomes strongly dependent on unknown 
properties of the surface states. In this paper we will not 
consider these close impurities, but discuss how details 
of the barrier may change the above result back to the 
e -(zi+z s )/a f R e f Qand thus allow the main contribution 
to come from pairs further from the barrier. 



Note also the extreme accuracy to which the positive 
and negative contributions to I r (z) cancel. Between two 
adjacent zeroes of Jo(k F r) the integral is of order 1 for 
small r (for r > a it gets damped by the exponential), 



-kpa 



io- 



yet the final integral is of order e kpZ 
if k F a = 100. 

A closer analysis of the integral over r shows that it 
accumulates a negative value of the order e~ z l a for small 
r, which is almost canceled by an equivalent positive con- 
tribution for very large r. The cancellation is found to be 
strongly dependent on the exact shape of the wave func- 
tion. Consequently, one may conclude that the cancella- 
tion can be lifted by specific properties of the tunneling 
amplitude, which could introduce an effective cut-off of 
the integration over r. In the following Sec. IIVI several 
models will be discussed, where the importance of hith- 
erto ignored details in the tunneling barrier will be made 
clear, and the importance of the assumed wave function 
for the localized states will be discussed. 

To complete the calculation of the conductance, we 
must now perform the summation over the Matsubara 
frequencies in the standard way, replacing 



r£/(iw n )= I 



de e 
/(e) tanh 



2T 



Under the assumptions we have made, this integration is 
independent of the details of the spatial integration, and 
will not be affected by the modifications introduced in 
later calculations. Integrating over the contour shown in 
Fig. |21 we get an expression for conductance depending 
on the discrete variable Q m . Making an analytical contin- 
uation of this function, and taking the limit as f£ m goes 
to zero, we finally get an expression for the conductance: 




FIG. 3: Integration contour. The values e = e a and e; — iQ m 
are shown by circles (o), the values e = ±A are shown by 
empty squares (□), while the poles of tanh(e/2T) are shown 
by filled squares. 



7r(27r) 4 a 2 e 2 |T |V 
2(k F a) 12 T 



G 



sin jk F p sl ) _ 2 Psi_ 



xF(z B )J r (zi)n(ei)n(e s )5(e s + ej) . 



(13) 
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Here n(e) = (e e / T + is the Fermi distribution. 

As shown in Ref. |jj it is important to include the effect 
of Coulomb interaction between the occupied sites which 
results in additional energy Uc = ^o^^P^i- However, ac- 
counting for this interaction does not change the strong 
cancellation. We will come back to the role of Coulomb 
interaction in Sec. where we discuss non-Ohmic con- 
ductance. In the next section we discuss how robust the 
cancellation is and how it is influenced by the properties 
of the tunneling barrier. 



IV. OHMIC CONDUCTANCE 

To understand how robust the cancellation shown in 
the previous section is, we will analyze several aspects of 
Ohmic transport through the interface. 



A. Importance of Impurity Wave Function 

The hydrogen-like wave function is a typical approxi- 
mation for shallow centers in semiconductors. However, 
the crossed Andreev transport can also take place in 
mesoscopic devices, e.g., between a bulk superconduc- 
tor and a pair of quantum dots. This is, in particular, 
the case for the previously suggested spin entangled 

If the above cancellation is correct also for that case, 
then the crossed Andreev transport would be hardly fea- 
sible since the dots would have to be located virtually at 
the interface. However, the tails of the wave functions 
of the electrons localized at quantum dots are far from 
being hydrogen-like. In general, they are dependent on 
the design of the quantum dots. In particular, for the 
lateral quantum dots designed by a properly engineered 
gate potential one can expect parabolic confinement. In 
this case the wave function tail is Gaussian rather than 
exponential. 

To check whether the above cancellation exist in this 
case we have repeated the calculations of Sec. IIIII but 
replacing # s (r) with <J G (r) = (2-Ka 2 )- 1 / 2 e-\ r - r ^ / 2a2 . 
As a result, the contribution of a given pair decays as 
e -(z s +zi)/a ^ anc j Qne re t urns to the estimates of Ref. 0. 
Thus the design of quantum dots chosen for the spin en- 
tangler can be crucial for it's potential success. 



B. Importance of Barrier Properties 

In Sec. lIIII we assumed that the tunneling amplitude is 
local and coordinate-independent, T(r,r') = Tq5(y — r'). 
It means that during tunneling an electron can transfer 
its momentum to some disorder-induced scatterers, and 
the tunneling amplitude is independent of the incident 
angle. To study the role of this simplification, we will 
proceed as follows. First, we consider the case of an ideal 
barrier for which the tunneling amplitude depends only 



on the incident angle. We will show that such depen- 
dence does not remove the cancellation, and the decay 
cx e ~( z i+ 2 s) fc F persists. Then we will consider the case 
of a barrier with fluctuating strength. We will find that 
fluctuations of a scale < a can strongly facilitate trans- 
port restoring the e~( zl+ZB >' a dependence. 

To make these consideration more specific let us as- 
sume that the effective barrier thickness, d, fluctuates 
along the interface. Then the tunneling amplitude is non- 
local, and the tunneling Hamiltonian acquires the general 
form of Eq. @. Consequently, the coordinates of the 
Green's functions for HI- and SC-side are different, and 
the proper diagram has the form of Fig. 0] rather than 
that of Fig. n For the following it is convenient to nor- 




FIG. 4: Left panel: The diagram of Fig. for the case of 
nonlocal tunneling transparency. Right panel: More detailed 
sketch of the coordinates. 

malizc the tunneling amplitude to T , 

T(r,r') = T /(r,r'). 

Then the spatial integral can be written as [cf. with 
Eq. ©] 

x y^F(r - n,u n )F^ (t -r 2 ,w„) 

xG(t', t' , uj n - Q m )G(T , 2 ,T , 1 ,u n ) . (14) 

Following the previous calculations, this can be split into 
separate, identical integrations for each impurity. To sep- 
arate the roles of barrier thickness fluctuations and de- 
pendence of the incident angle let us express the tunnel- 
ing amplitude as 

f(T,T')= 9 (T)h(T'-T). (15) 

where g(r) describes spatial variations in the barrier, 
while h(r' — r) accounts for the dependence on the in- 
cident angle. Here both vectors r and r' belong to the 
interface plane. The function h(r' — r) can be assumed 
to depend only on |r' — r|. 
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a. Smooth Barrier: Let us start with the case when 
g(r) = g = constant. The basic spatial integration is 

I(r,r ,u n ) = j ' d 2 r"G{\r + v"\,r' ,u n )h{r"). (16) 

Here we have taken into account that the Green's 
function G(r, ri,w n ) depends only on V z 2 + r 2 and 
sj z 2 + r 2 . Let us now assume that h(r") decays much 
more rapidly than G. That allows us to expand the in- 
tegrand in powers of x" and y" keeping only the second 
order (the first order term vanishes on integration), 



G(|r- 



G 



d 2 G 



d 2 G 



2 dx 2 



2 dy 2 



with G = G(r,r' ,tu n ). Now let us consider the simplest 
case of a rectangular barrier for which the function h(r) 
can be modeled as h(r) = d~ 2 e~ B< - r / d ^ where B > 1 is 
some dimensionless constant. This model follows from an 
assumption that the barrier is uniform along the surface 
and rectangular. Then the tunneling exponent can be 
written as 



dJ2m[U - £7(1 - k 2 /k 2 )} « -A - B{r/df 



where A = dyj2m{U - E), B = \ E/2(U-E), U is the 
barrier height, while E is the electron energy. Using this 
model and explicit expression |(HJ| for the Green's function 
we obtain 



c(r) 



B 



I(r, r Q ,uj n ) 
d 2 I 

1 



G(r,r' ,oj n 

,.2 



c(r) . 
2a 



ar 



Vr 



+ z z Vr 2 + z 2 



The integral I r , Eq. 1|12[) . now changes to 



Ir{z) 



rdrJo(hFT)e 



-y/r 2 -\-z' 2 / a 



c(r). 



At d <C a the correction induced by the angular depen- 
dences of the transmission is small, but, in principle, it 
could be sufficient to lift the cancellations that give the 
e _fcF2 -dependence. However, a combination of analytical 
and numerical analysis shows that the additional terms 
also lead to the e _fcFZ -dependence. We therefore con- 
clude that the dependence of the tunneling transparency 
of a uniform barrier on the incident angle still leads to 
the e _fcFZ -decay of the crossed Andreev transport. 

b. Inhomogeneous barrier: We now turn to the sit- 
uation were we have fluctuations in the barrier strength 
that facilitate tunneling through the places where the 
barrier is thin. Let us assume that the typical size of 
these regions, A is much less than the localization length, 
a, but larger than kp 1 , 

kp 1 < A < a. 

Now we cannot assume g(r) to be constant over the re- 
gion spread by the impurity potential. For simplicity, we 



will in this case write h(r") — S(r"), as we have pre- 
viously shown that these corrections do not change the 
principal behaviour of the transport. The shape of g(r) 
is dependent on the roughness, and on the relative po- 
sitions of the barrier minimum and the impurity center. 
For a barrier with a parabolic minimum one can show 
that g{r) has a Gaussian shape, 



9( r ) 



2 / 2 

/a 



The general analysis of this situation is complicated, as 
the impurity center may not coincide with the center of 
the barrier minimum. Two simplified cases arc still suf- 
ficient to shed light on the situation. 

Let us for simplicity start with the case when the mini- 
mum in the barrier strength coincides with the projection 
of the impurity center on the interface. In this case the 
integral for I r (z) similar to Eq. (|12|l can be analyzed in 
detail. It turns out that with increasing z it crosses over 
from e~ kpZ to e~ z ' a at some z* which depends on a. 
The quantity z* decreases with decrease of a, z* = a at 
1/(7 = 0.15. Thus the barrier inhomogencity (modeled 
by small a) facilitates the transport by eliminating the 
cancellation. Now we can relax the previous assumption 




FIG. 5: The dependence I r (z) for fcfffl = 20 and \ja = 
0.01,0.1,0.15,0.2. 

that the barrier strength minimum occurs exactly at the 
projection point of an impurity center and consider the 
situation where the minimum is off-center, but varies so 
fast that G(r) can be considered constant in comparison. 
In this case we can replace G(r) = G(r m j n ), r mm being 
the barrier minimum coordinate, and move the origin of 
the integration to r m i n . The integral will then again give 
a simple e~ 2 / a -dependence. Thus we see that the cut-off 
introduced by a clear minimum in the barrier is suffi- 
cient to change the ^-dependence of the conductance per 
impurity. 

If there are several minima within one single impurity, 
one could imagine that these could give new interference 
effects. However, integration over several minima cor- 
responds to simply summing up the contribution from 
the separate integrals. Each minimum will be coupled 
to a minimum on another impurity, giving a prefactor of 
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sin 2 p s ;. When several sines are added at each impurity, 
we get 

(sinp s z +sinp s /;/) 2 = sin 2 p s ; + sin 2 /v// +2sinp s ; sin/V!'- 

When averaging over several pairs, the sin 2 -terms will 
survive, while the cross-terms will average to zero. We 
therefore assume several minima within the range of each 
impurity to be equivalent to several separate pairs in the 
total averaging. Considering the probability of a pair 
accepting a Cooper-pair, the effect on one minimum of 
the pair already being occupied due to another minimum 
should be negligible. 

The characteristic localization length of the electron 
under the barrier, we call a. In order for these consider- 
ations to be relevant, the barrier thickness d must vary 
with several a over a length scale much shorter than a. 
Changes in barrier thickness that do not meet this condi- 
tion are better analyzed in terms of the following model. 

c. Barrier with a block-like disorder: The model dis- 
cussed above relies on a change in barrier thickness that 
is of the order of a. Because of the very fine cancella- 
tions in the integral I r , Eq. I|12(l . much smaller changes 
in barrier height can be important, provided they are on 
a length scale of the order of a. This can be demonstrated 
using a simple model based of the analysis of the integral 

I r {R,z) = / rdrJ (k F r)e~ Vr ' +z " /a 
Jo 

as a function of the cut-off R. Obviously, I r (z) = 
lirriR^oo I r (R, z). By splitting the integration over x into 
intervals divided by the zeroes of the Bessel function, 
we get successive contributions of alternating signs and 
close absolute values. The result approaches the e~ kpZ 
behaviour seen before when X goes to infinity. The ab- 
solute value of the sum, S(n, z) of an even number of the 
intervals defined as 

2 " rR m+ 2 

S{n,z) = rdrJ {k F r)e- Vr " +zA/a 
m =0 *'- R m 

= f 2 " +2 rdrJ (k F r)e- Vr ' +z ^ a 
Jo 

where R m /a is the m'th zero of the Bessel function, 
Rq = 0, will therefore have a maximum for some n, as 
shown in Fig. EI This maximal value of S(n, z) this max- 
imum varies with z as e _z / a , and the corresponding cut- 
off radius, R* = i?n max , is a slowly varying function of 
the ration z/a, corresponding to 20 to 100 zeroes of the 
Bessel function. 

Based on this property we construct a simple model of 
a barrier with a block-like disorder assuming 

g(r) = l + r]Q(R* -r) (17) 

where O(r) is the Heavisidc step-function. Since this con- 
tribution only slowly decays with z, even a small barrier 
variation, i], can give significant contributions. Since the 
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FIG. 6: S(n), the integral cut-off after 2n zeroes of the Bessel 
function, k F a = 100, z/a = 1 



first term in Eq. (|17fl leads to a decay oc e~ kpZ while 
the second contribution decays oc e~ z l a we only need 
rj > e -(fci ? - 1 / ) z . Figure [7| shows the optimal cut-off ra- 
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FIG. 7: The functions n m ax(z) (top line) and I r (R*,z) (mid- 
dle line) for k F a — 100. The bottom line shows the function 

e~ z/a . 

dius, R* , as a function of z, as well as I r {R* , z) compared 
to a graph showing e~ z / a . 

Obviously, small variations in inhomogeneity size, R, 
will give large variations in the result, this can be reme- 
died by smearing the cut-off over a period or two of the 
Bessel function. If we assume the center of the barrier 
reduction to be slightly displaced from the impurity cen- 
ter, the integration over angles can probably be found to 
be some such smeared step function. At some point, the 
smearing will be such that the cancellation is no longer 
lifted, and the barrier minimum gives insignificant con- 
tribution. 

Thus we see that a small variation in the barrier thick- 
ness, as long as it is at the proper length scale and cen- 
tered on the impurity center, can drastically reduce the 
interface resistance of the barrier. 
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V. NON-OHMIC CONDUCTANCE 

Having established that there exists a range of appli- 
cability of the underlying model we now address the non- 
Ohmic behavior of the interface conductance. To figure 
out the nonlinear properties one has to compare eV with 
other relevant energy scales: temperature T, inter-site 
Coulomb repulsion energy Uc, and the energy splitting 
e s — ei of a pair. In addition, we have to consider the 
width of the <5-function in energy that selects which pairs 
may contribute. This width can be estimated as a typical 
inverse life time of an electron at a localized site forming 
the pair. We always assume eV <C A excluding in this 
way the possibility of single-electron transport. 

The net current over the barrier can be seen as the dif- 
ference between current from the superconductor to the 
insulator, Is-*i and current from insulator to supercon- 
ductor Ii-ts- These currents are in turn determined by 
the transition probability and the occupation probabil- 
ities for the involved states. Following Ref. [ll| we will 
assume all matrix elements to be energy- independent. 
Thus the only variations in the transition rate is due to 
the occupation probabilities. We must also remember the 
eV-dependence of the (^-function, describing the conser- 
vation of energy, which regulates which pairs contribute 
to the transport. 

Since the superconducting condensate has a macro- 
scopic number of states at the level of the chemical poten- 
tial, the current Is—ti is only dependent on the probabil- 
ity of finding an empty pair in the insulator, while If_>s 
requires an occupied pair. In both cases the relevant pair 
will have to satisfy energy conservation. 

To keep track of realistic situations we assume that 
the entire voltage drop occurs at the barrier, but al- 
low for a small portion of the insulator near the inter- 
face to be filled up or emptied by electrons due to the 
voltage drop over the barrier. This region models the 
depletion zone of a semiconductor heterojunction. For 
impurities outside the depiction zone, the Fermi level is 
assumed to be fixed relative to the impurity energy lev- 
els. In this case the (5-function in energy must be chosen 
as 5(e s + ei + Uc — 2eV), and the occupancy numbers are 
given as before, as n(e). Very close to the barrier, inside 
the depletion layer, we instead use a picture where we 
keep the impurity energy levels constant relative to the 
superconductor condensate, but adjust the Fermi level to 
get n(e, eV) — n(e + eV). 

eV > eV = eV < eV > eV = eV < 







6/ 










€ s ° 




— e,* 


£/ 




£° 


£(° 











FIG. 8: Electrochemical potential in HI (left) and SC (right) 
outside (left panel) and inside (right panel) the depletion zone 



In both cases we have to consider the Coulomb energy 
Uc, and the simplest way of accommodating it is by de- 
scribing each pair as a four-level system corresponding to 
the 4 following configurations: (i) both sites are empty, 
(ii, iii) one site is occupied, and (iv) both sites are occu- 
pied. The configurations are shown in Fig. El For sim- 



FIG. 9: Energy levels for 4 configurations of the pair: e 1 — 0, 

£10 = 6s, £01 = 61, £11 = £s + 6; + Uc- 

plicity we disregard Coulomb interaction with charges 
outside the pair. We can then write a partition func- 
tion for the four-level system, and use this to find the 
probabilities of a pair being empty - allowing a Cooper 
pair to fill it - or if filled, allowing one Cooper pair to 
be created. Using the energies defined in the figure, the 
partition function can be written as 

Z =1 + e -( £ -'-^)/ T + e -(ei-»)/T + e -(e s + £i -2 M +(7 c )/T 

and the probabilities of the different configurations in 
similar notation: 

Poo = Z~\ P W = Z-'e-^~M T , 

P l = Z-V^-^/r Pu = Z -l e -(e 3+£l -2 M+ l/ c )/T 

Here fi is the chemical potential in the superconductor. 

The current through the interface is proportional to 
the difference 

Poo - Pu = Z- 1 (e-(*s+ei-^+Uc)/T _ A (18) 

With the inclusion of the Coulomb energy this can no 
longer be factorized into separate occupation probabili- 
ties of the two impurities of the pair, but if we set Uc = 0, 
it can be seen that the difference simplifies to the former 
results. 

If we consider the pairs inside the depletion zone, we 
use fi = —eV. Taking into account the energy conserva- 
tion law requiring e s + e; + Uc = one can express the 
difference Poo — Pu as 

I _ e 2eV/k B T 

1 + e (-e s +eV)/fe B T + e (e B+Uc +eV)/k B T + e 2eV/k B T 1 

In this case, the energy conservation law is independent 
of eV, so the transitions are suppressed until the voltage 
reaches |e s | + Uc, then it rapidly rises, before Poo — P\i 
saturates at unity. For the pairs outside the depletion 
zone, we can write fx = and the expression for Poo — Pu 
is 

e 2ey / T - 1 

1 + e - e */ T + e (^+2eV+U c )/T _|_ e 2eV/T ' ' 

While these expressions are somewhat similar, the main 
difference lies in the el^-dependence of the conservation 
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law which means that as the voltage changes, the choice 
of pairs satisfying energy conservation will change, so a 
single pair will pass into and out of the allowed range 
instead of saturating. 

In the following we assume the Fermi level in a semi- 
conductor to located inside the impurity band, which is 
sufficiently larger than A 3> eV. Therefore when per- 
forming the summation over all pairs, we choose a uni- 
form distribution of e s . In this case, both expressions will 
give exactly the same results, although the physics be- 
hind them are slightly different. Neglecting the Coulomb 
interaction between the components of the pair, we sim- 
ply get a linear relation in both cases. Thus Ohmic be- 
havior would persist even though eV can exceed the tem- 
perature T. For pairs inside the depletion zone, this is 
due to inclusion and saturation of more pairs as eV be- 
comes larger than the energy splitting of the pairs. If a 
given pair has started contributing, one channel has been 
opened, it will not terminate with increasing eV . Out- 
side the depletion zone, pairs will only contribute for the 
width of the (5-function, so as voltage rises, other pairs 
will take over the transport, but the number of pairs that 
can contribute will rise linearly with the voltage. 

A typical / — T^-curve is shown in Fig. 1101 This curve 
is calculated under a simplifying assumption that the 
Coulomb correlation energy Uc is kept constant of the 
same order as the temperature T since the contribution 
of the si pair is cx pJ; 2 e~ Ps! ^ and cut-off at small p by 
the requirements Uc ~ For this pairs the Coulomb 
interaction is essentially screened and does not block two- 
electron tunneling. This assumption significantly simpli- 
fied the calculation comparing with averaging over all 
correlation energies, but does not change the conclusion. 
As we see, the transport is suppressed at eV < l^c, while 
for eV > Uc transport will return to a Ohmic behavior, 
as shown in Fig. 1101 Thus we predict a rather unusual 
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Coulomb correlation between the occupation numbers of 
the pairs responsible for crossed Andreev reflections. 
VI. CONCLUSIONS 

In conclusion, we have shown that the pair tunneling 
through a barrier at the interface between a supercon- 
ductor and hopping insulator is extremely sensitive to the 
properties of the tunneling barrier. This sensitivity is due 
to rapid oscillations (at scale ~ kp 1 ) of the electron wave 
functions in a superconductor comparing to the charac- 
teristic scale, a, of variation of the localized wave function 
in a hopping insulator. These oscillations cause dramatic 
cancellations in the tunneling probability if the barrier is 
uniform. However, this cancellation of the interference 
origin can be suppressed if the barrier is inhomogeneous, 
as it was demonstrated for different models of a barrier. 
In particular, the fluctuations in the barrier strength of 
the scale A falling within the window 1/kp 'ti X < a sup- 
press the cancellations and restore the transport even if 
their relative amplitude rj is very small. For a barrier 
with a block-like disorder with the scale ~ a we obtained 
an estimate for suppression of the oscillations r\ > e~ kFa 
for the pairs located at the distance ~ a from the in- 
terface. This happens if the beneficial barrier fluctua- 
tions must coincide with particularly positioned impuri- 
ties with the right energies, means that the number of 
impurities contributing to conductance will be relatively 
small and the relative importance of "successful" pairs 
will increase. Consequently, one can expect pronounced 
mesoscopic - sample specific and reproducible - fluctua- 
tions in both Ohmic and non-Ohmic conductance. Such 
fluctuations will be especially pronounced when the bar- 
rier contains large transparency fluctuations (punctures). 
One can expect that the fluctuations will have different 
behavior depending on the location of the relevant pairs 
with respect to the position of the depletion zone near 
the interface. We plan to study mesoscopic fluctuations 
of the Andreev transport between a superconductor and 
a hopping insulator in more detail as a separate project. 

Another specific feature of the transport is sensitiv- 
ity of the non-Ohmic transport to Coulomb correlation 
in the occupation numbers of the relevant pairs. This 
correlation leads to non-Ohmic behavior at low voltages, 
eV < Uc, T, while at higher voltages the transport turns 
out to be Ohmic. This re-entrant behavior is a hallmark 
of the Coulomb correlation. 
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FIG. 10: Current-voltage curve for Uc /T = 0, 1, and 2. 

situation when the conductance evolves for a low-field 
Ohmic to a high-field Ohmic regime through an interme- 
diate non-Ohmic one. Such behavior is a hallmark of the 
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